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Abstract 



'^ ■ We study a contact process with creation at first- and second-neighbor sites and inhibition at 

first neighbors, in the form of an annihilation rate that increases with the number of occupied first 
__i , neighbors. Mean-field theory predicts three phases: inactive (absorbing), active symmetric, and 

Q ' active asymmetric, the latter exhibiting distinct sublattice densities on a bipartite lattice. These 

cj : 

phases are separated by continuous transitions; the phase diagram is reentrant. Monte Carlo simu- 
lations in two dimensions verify these predictions qualitatively, except for a first-neighbor creation 
rate of zero. (In the latter case one of the phase transitions is discontinuous.) Our numerical 



in 

(^ • results confirm that the symmetric-asymmetric transition belongs to the Ising universality class. 
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and that the active-absorbing transition belongs to the directed percolation class, as expected from 
symmetry considerations. 
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I. INTRODUCTION 

An absorbing-state phase transition is a far from equilibrium transition between an active 
and an inactive phase. When a control parameter such as a creation or annihilation rate is 
varied, the system undergoes a phase transition from a fluctuating state to an absorbing one. 
Such transitions are currently of great interest, arising in a wide variety of problems, such 
as population dynamics, heterogeneous catalysis, interface growth, and epidemic spreading 

5|. Interest on this class of phase transitions has been stimulated by recent experimental 



realizations 



BE 



As equilibrium statistical mechanics, absorbing-state phase transitions can be classified 
into universality classes. However, a complete classification of their critical behavior is still 
lacking; much numerical and theoretical work is focused on identifying the factors that 
determine the universality classes of these transitions. 

The absorbing-state universality class associated with directed percolation (DP) has 
proven to be particularly robust. DP-like behavior appears to be generic for absorbing- 
state transitions in models with short-range interactions and lacking a conserved density or 



a. 



symmetry beyond translational invariance [8[. By contrast, models possessing two absorb 



la 



ing states linked by particle-hole symmetry belong to the voter model universality class [9|. 
There are also models which, while free of absorbing states, cannot achieve thermal equilib- 
rium because their transition rates violate the detailed balance principle. An example is the 
majority vote model. In its ordered phase a Z2 symmetry is spontaneously broken, leading 
to Ising-like behavior in two or more dimensions. 

In this work we examine whether a spatially structured population model that suffers a 
phase transition to a single absorbing state can also exhibit a broken-symmetry phase. Our 
candidate for such a model is a modified contact process (CP) with suppression of activity at 
the nearest neighbors of active sites. We observe unequal sublattice populations for suitable 
choices of the control parameters. 

The balance of this paper is organized as follows. In the next section we define the model 
and analyze its mean-field theory. In Sec. Ill we present our simulation results; Sec. IV is 
devoted to discussion and conclusions. 



II. MODEL AND MEAN-FIELD THEORY 



The CP [lO| is a stochastic interacting particle system defined on a lattice, with each site 
i either occupied by a particle [(Tilt) = 1], or vacant [crj(t) = 0]. Transitions from ctj = 1 to 
(Tj = occur at a rate of unity, independent of the neighboring sites. The reverse transition, 
at a vacant site, is only possible if at least one of its neighbors is occupied: the transition 
from (Tj = to cr^ = 1 occurs at rate Ar, where r is the fraction of nearest neighbors of site 
i that are occupied. In the absence of spontaneous creation of particles, the state cr^ = 



for all i is absorbing. It can be shown 10|] that for a certain critical value Ac the system 



undergoes a phase transition between the active and the absorbing state. Although no exact 
results are available, the CP has been studied intensively via series expansion and Monte 
Carlo simulation, and its critical properties are known to high precision [l|-l5|. 

In order to generate distinct sublattice occupations, we modify the basic contact process 
in two respects: 

i) In addition to creation at first neighbors, at rate Ai, we allow creation at second 
neighbors, at rate A2. For bi-partite lattices such as the square or simple cubic lattice, Ai is 
the rate of creation in the opposite sublattice, while A2 is the rate in the same sublattice as 
the replicating particle. Unequal sublattice occupancies are thus favored for A2 > Ai. 

ii) We include a nearest-neighbor inhibition effect in the annihilation process. In addition 
to the intrinsic annihilation rate of unity, there is a contribution of finf, where ui is the 
number of occupied first neighbors. Thus if the occupation fraction pA of sublattice A is 
much larger than that of sublattice B, any particles created in sublattice B will die out 
quickly, stabilizing the unequal sublattice occupancies. 

The one-site mean-field theory (MET) for this model on a lattice of coordination number 
q is given by the coupled equations, 

^ = -(1 + pq^pDpA + (Aips + X2Pa){1 - Pa) (1) 

and 

-(1 + pq^p\)pB + (AiPA + A2Pb)(1 - Pb) (2) 



dt 

which are seen to be symmetric under pA ^ Pb- Defining p = Pa + Pb and (p = Pa — Pb, 
these equations may be recast as. 



dp 
~(E 



(A-I)P-^P^ 



y0^ - ^pg^(/ - <P')P 



and 



dt 



A 



1 - A2P - ^/ig'(p' - 0') 



(3) 



(4) 



where A = Ai + A2 and A = A2 — Ai. 

From Eq. ([3]) it is evident that at this level of approximation, the extinct ion- survival 
transition occurs at A = 1, as one would expect. This transition leads to the symmetric 
active phase, characterized by p > and = 0. This phase is linearly stable for small A. As 
A2 is increased, with Ai constant, the symmetric active phase can become unstable, leading 
to a phase with |0| > 0. On further increasing A2, however, p can increase to the point 
that the symmetric active phase is again stable. (Perfect sublattice ordering corresponds to 
p = 0, which is only possible for p < 1; total densities larger than unity tend to suppress 
sublattice ordering.) Thus the asymmetric active phase should exist for some intermediate 
range of A 2 values. 

In the symmetric active phase, the stationary activity density is. 



2k L 



v/(A/2)2 + 4«:(A-l) - (A/2) 



(5) 



with hi = fJ^q'^/4:. From Eq. (^, this solution is stable with respect to one with 7^ 0, when 



a^ = A — 1 — A2P + Kp^ < 0. 



(6) 



The transition to the asymmetric phase occurs when 



0, which implies 



(2A2 + A)2(A - 1) = 2(A2 - 1) [8fi:(A2 - 1) + (2A2 + A)A] 



(7) 



Eq. (jl]) may be written as d(j)/dt = a^(f) — K,(j)^, showing the expected symmetry under — )■ 
—0, which in turn suggests that the transition between symmetric- and asymmetric-active 
phases is Ising-like. The two transitions (active/absorbing and symmetric/asymmetric) co- 
incide at the point Ai = 0, A2 = 1. Figure [1] shows the phase diagram for p = 2, obtained 
using Eq. (^^. Of particular note are the observations that (1) the asymmetric-active phase 
is observed only for intermediate values of A2, i.e., the phase diagram is reentrant, and (2) 




FIG. 1: (Color online) Phase diagram in the Ai — A2 plane for fi = 2, showing absorbing (ABS), active- 
symmetric (AS) and active asymmetric (AA) phases. Lines: MFT; symbols: simulation. 



the asymmetric-active phase does not exist for Ai above a certain value, A|(/i). For /i = 2, 
for example, A^ = 7.1443. Also shown in Fig. [T]are simulation results (details of our simula- 
tion method are discussed in the following section). MFT and simulation are in qualitative 
agreement, but as is usually the case, MFT overestimates the regions in parameter space 
corresponding to the active and ordered phases. In particular, MFT overestimates AJ by 
about an order magnitude. 

The line Ai = is special, since the subspaces with p^ > and pb = (and vice- versa), 
are also absorbing. MFT yields the following phases: absorbing for A2 < 1; asymmetric- 
active-I for 1 < A2 < A'^(/i); asymmetric-active-II for A^(/i) < A2 < \^^{fi); and symmetric- 
active for A2 > A^^. The difference between phases asymmetric-active-I and II is that in the 
former, p = (i.e., only one of the sublattices is populated) while in the latter p > > 0. 
The transitions between these phases are continuous. (Naturally, the asymmetric-active-II 
and symmetric-active phases can only be reached from initial conditions with both sublattices 
populated.) For p = 2, for example, we find A^ = 30.97 and A^^ = 41.65. 



We also studied a version in which the suppression (enhanced annihilation rate) is a linear 
function of the number of occupied first neighbors, i.e., the mean- field equations 

—A = _(i + fiqpB)pA + (AiPb + A2Pa)(1 - Pa), (8) 

and similarly for Pa ^ Pb- The phase diagram is qualitatively similar to that found above, 
but the region occupied by the asymmetric-active phase is much reduced. For /i = 2, for 
example, sublattice ordering is only found for Ai < 0.6805. Since mean-field theory offers at 
best a qualitative description of the model, we turn, in the following section, to simulation. 

III. SIMULATIONS 

We performed extensive Monte Carlo simulations of the model on square lattices of linear 
size L = 20,40, ..., 320 sites, with periodic boundaries. The simulation algorithm is as 
follows. First, a site is selected at random. If the site is occupied, it creates a particle at 
one of its first-neighbors with a probability pi = Ai/(1 + Ai + A2 + P'nf), or at one of its 
second-neighbors with a probability p2 = ^2/i^ + Ai + A2 + pnf). With a complementary 
probability 1 — {pi + P2) the site is vacated. (In order to improve efficiency the sites are 
chosen from a list which contains the currently Nocc occupied sites; we increment the time 
by At = 1/Nocc after each event). For simulations in the subcritical and critical absorbing 



regime, we employ the quasi-stationary (QS) simulation method detailed in [ll|, to further 
improve efficiency. 

In a series of studies using /i = 2.0, we verify that for low values of Ai and A2 the 
system becomes trapped in the absorbing (inactive) state. When we increase Ai and/or A2 
the system undergoes a phase transition from the inactive to the active symmetric (AS) 
phase. For example, for Ai = 0.2, the absorbing transition occurs at A2 = 1.5620(5). (Note 
that at this point Ai -|- A2 is somewhat greater than the critical value for the basic process. 



Ac = 1.6488(1), on the square lattice [12[.) At the critical point, the quasistaionary order 
parameter is expect to decay as a power law, p ~ L~^/'^-^. In Fig. Illlt we show that 
our data follow a power law, with /3/z/^ = 0.79(1), in good agreement with the DP value 
(5 /u±_ = 0.797(3) [13|. The lifetime of the QS state r ~ L'^W''^^ at criticality, where we find 
^\\/^± = 1.75(2), also shown in Fig. 2. This value again is in very good agreement with the 
best known DP value lyw/iy^ = 1.7674(6). The moment ratio m = (p^)/(p)^ (see Fig. [3]) is 

6 



analogous to Binder's fourth cumulant [ij] at an equilibrium critical point, in the sense that 
m assumes a universal value rric at the critical point. We find = 1.324(5), in comparison to 
the known DP value m = 1.3264(5) [l5(]. These results permit us to state that, as expected, 
the absorbing transition belongs to the universality class of directed percolation. 




InL 



FIG. 2: (Color online) Scaled critical QS density of active sites In p (bottom) and scaled lifetime of the QS 
state In t (top), versus InL. Parameters: /.t = 2.0, Ai = 0.2 and A2 = 1.5620 




FIG. 3: (Color online) Quasistationary moment ratio m versus A2, for fi = 2.0 and Ai — 0.2. System sizes: 
L = 20, 40, 80, 160 (in order of increasing steepness). 



In the active phase, we find that as we increase A2 at fixed Ai, there is a transition to 
a state with unequal particle densities in the two sublattices (</> 7^ 0), that is, to the active 
asymmetric (AA) phase. A typical configuration in the AA phase is shown in Fig. HI Our 
observation of a bimodal probability distribution for the order parameter (p (see Fig. IIIip 
confirms the existence of the AA phase. Although MFT predicts a symmetry-breaking 
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FIG. 4: (Color online) A typical configuration in the active asymmetric phase; particles are represented by 
black sites. System size: L = 50. Parameters: /i = 2.0, Ai = 0.2, and A2 = 6.0. 
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FIG. 5: Probability distribution of the order parameter cj). System size: L = 80. Parameters: // = 2.0, 
Ai = 0.2. A2 = 2.0(dotted curve), 10.0 (dashed) and 40.0 (solid). 



transition in any number of dimensions, we do not expect such a phase transition in one 
dimension; simulations on a ring confirm this. 

In order to classify the AA-AS transition, we perform a finite-size scaling study. First we 



investigate the behavior of the reduced Binder cumulant [IJ] , given by 



Q, 



Ua = 1 



2\2' 



(9) 



The intersection points of the cumulants for successive pairs of sizes depend rather weakly 
on the sizes, providing a good estimate for the critical point. The value of the cumulant at 
the critical point approaches a universal value as L — t- oo. In Fig. [6] we show the results for 
the case /i = 2.0 and Ai = 0.2. The curves for different sizes intersect at A2 = 3.940(5), and 
again at A2 = 31.92(6) when L — t- 00. 

Extrapolating to infinite size, we estimate f/4c = 0.605(10) at the transitions, consistent 
with the Ui^c = 0.61069... 16| for the two-dimensional Ising model with fully periodic bound- 
ary conditions. For values of A2 between the transitions points, the cumulant approaches 
2/3, as expected in an ordered phase that breaks a Z2 (i.e., up-down) symmetry. 

Further evidence for Ising-like behavior is furnished by the order-parameter variance. At 
the critical point, var(0) = (0^) — (0)^, is expected to scale so: L^var(0) oc L'^''^; this is 
confirmed in Fig. [71 A fit to the data yields the exponent ratio '-f/u = 1.76(5), again in good 
agreement with the known value of 7/1/ = 7/4 for the Ising model in two dimensions. 




FIG. 6: (Color online) Binder cumulant versus A2, for fi — 2.0 and Ai — 0.2. System sizes: L — 20, 40, 80, 160 
(in order of increasing steepness). 



For the Ai = 0, simulations show no evidence of the asymmetric-active-II phase predicted 
by MFT. Simulations reveal only two phase transitions along this line. The first, from 




FIG. 7: (Color online) Scaled order parameter variance InL^ var0 versus InL at the AS-AA transition. 
System sizes: L = 20, 40, ..., 320. Parameters: /i — 2.0, Ai = 0.2. A2 = 3.94. The slope of the regression line 
is 1.76. 

the absorbing to the asymmetric active phase, occurs, as expected, at the critical value 
of the basic CP on the square lattice. We find A2,c = 1.6488(3), and confirm that this 
transition belongs to the DP class (see Fig. [HI where we find m = 1.320(8), (S/u^ = 0.81(1) 
and i^\\/i^i_ = 1.74(3), consistent with the values for DP in 2+1 dimensions, cited above). 
The second transition, between the active-asymmetric and the active symmetric phases, is 
discontinuous. Figure [9] illustrates such transitions for /i = 0.2 and /i = 2. These results 
signal a qualitative failure of mean-field theory along the line Ai = 0. 




FIG. 8: (Color online) Quasistationary moment ratio m versus A2, for fi = 2.0 and Ai ~ 0. System sizes: 
L = 20,40,80, 160 (in order of increasing steepness). Inset: Scaled critical QS density of active sites In p 
(bottom) and scaled lifetime of the QS state In r (top), versus InL. Parameters: p. = 2.0, Ai = and 
A2 = 1.6488 
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FIG. 9: (Color online) Order parameters p (solid) and (\) (dashed) versus A2, for Ai 
/i = 0.2 and /i — 2.0 from top to bottom. System size: L=80. 
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IV. CONCLUSIONS 

A contact process with creation at both nearest and next-nearest neighbors, and suppres- 
sion (in the form of an increased annihilation rate) at neighbors of active sites is shown to 
exhibit a broken-symmetry phase with sublattice ordering. The existence of a continuous 
symmetry-breaking transition is predicted by mean-field theory and confirmed in simula- 
tions on the square lattice. The asymmetric active phase exists for a restricted range of the 
second-neighbor creation rate, A2, when that at first neighbors, Ai, is sufficiently small. For 
a given value of the suppression parameter /i, there is a certain value, X^, above which no 
sublattice order occurs. We note that the possibility of breaking Z2 symmetry depends on 
the fact that creation rates Ai and A2 apply to different sublattices. We would not expect 
to observe this symmetry breaking if, for example, there were a creation rate Ai for both 
nearest- and second-neighbors, and a different creation rate at third-neighbor sites. 

For < Ai < A*, the sequence of phases observed on increasing A2 from small values, while 
holding /i and Ai fixed, is: absorbing, symmetric-active, asymmetric-active, and symmetric- 
active. The first transition belongs to the universality class of directed percolation, while 
the latter two are Ising-like. For the special case of Ai = 0, the first transition is from 
absorbing to asymmetric-active and the second is discontinuous. For both values of n studied 
we observe continuous transitions for Ai > 0. While increasing Ai makes the transition 
smoother, we can not exclude the possibility of a discontinuous transition when Ai > 0, 
for larger values of /i. A complete characterization of the reentrant transition for yU > 2.0, 
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including the location of the possible tricritical points separating the line of second order 
transitions from that of first order transitions is left for future work. 

The absorbing- active transition at Ai = marks the meeting of two lines of critical 
points, one of DP transitions, the other of Ising transitions. The meeting of two such lines 



has been studied in the context of a nonequilibrium Potts model by Droz at al. [17[, and 
of the generalized voter model (GVM) by Al Hammal et al. [18|. The former study also 
finds an absorbing state in the DP class, while in the latter, after the Ising and DP lines 
join, the absorbing transition belongs to the GVM class. We note that in these studies the 
models have but two (symmetric) absorbing states (i.e., all plus and all minus), whereas 
for Al = 0, our model has two (symmetric) absorbing subspaces (activity restricted to one 
sublattice) and a third, completely inactive absorbing state. In terms of long-time behavior, 
the absorbing-active phase transition takes the system from the fully inactive state to one 
with activity restricted to one sublattice. Staring from all sites active, at the critical point 
Al = 0, A2 = Ac,cP; fluctuations drive local regions into states with activity in only one 
sublattice, as well as fully inactive regions. We expect that these domains then coarsen, so 
that eventually there is activity only in one sublattice; from then on the dynamics is that of 
the basic contact process. Thus it is not surprising that we observe DP-like static behavior. 
It is nevertheless possible that the short-time critical behavior includes a regime dominated 
by domain dynamics as in the GVM. It also seems possible that the discontinuous AA-AS 
transition observed for Ai = belongs to the voter model class. We intend to explore these 
questions in future work. 

Possible extensions of this work include precise determination of static and dynamic crit- 
ical properties at the symmetry-breaking transition, studies of the effect of diffusion and/or 
anisotropy (which may generate modulated phases) and studies of a model in continuous 
space. The latter appears to be of particular interest in the context of patterns arising in 
ecosystems or bacterial colonies. We note that the suppression of activity at the nearest 
neighbors of active sites resembles biological lateral inhibition, known to be important in 



the visual system of many animals [19[. Thus extensions of the present study might be 
useful in the study of patterned activity on the retina. Finally, if the inhibition rate were a 
decreasing function (linear or quadratic) of the number of occupied nearest neighbors (i.e., 
yU < 0), one should observe symbiotic coexistence of activity in the two sublattices. The 
latter is reminiscent of the Janzen-Connel hypothesis for the maintenance of tree species 
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biodiversity in tropical rainforests 20 1. 
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